In [1]:
%matplotlib inline
#%matplotlib notebook
%load_ext version_information
%load_ext autoreload
In [51]:
import datetime
import os
import sys
import warnings
warnings.simplefilter("ignore")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import spacepy.datamodel as dm
import spacepy.toolbox as tb
import tqdm
import pymc3 as mc3
import seaborn as sns
sns.set()
%version_information matplotlib, numpy, pandas
Out[51]:
In [76]:
plt.rcParams['figure.figsize'] = [10, 6]
plt.rcParams['savefig.dpi'] = plt.rcParams['figure.dpi'] # 72
%config InlineBackend.figure_format = 'retina'
In [81]:
# Lets make to Normal distributions and add them together
locs = np.asarray([3,10])
sds = np.asarray([2,2])
D1 = np.random.normal(locs[0], sds[0], size=100)
D2 = np.random.normal(locs[1], sds[1], size=80)
In [82]:
# bins = np.linspace(-24,24, 100)
sns.distplot(D1, kde=True, label='D1')
sns.distplot(D2, kde=True, label='D2')
# plt.hist(D1, bins=bins, histtype='step', lw=3, label='D1')
# plt.hist(D2, bins=bins, histtype='step', lw=3, label='D2')
D = np.append(D1,D2)
sns.distplot(D, kde=True, label='D1+D2')
plt.legend()
Out[82]:
In [113]:
D.shape
Out[113]:
In [114]:
SEED = 20161210
with mc3.Model() as model:
w = mc3.Dirichlet('w', np.ones_like(locs))
mu = mc3.Uniform('mu', -10., 20., shape=locs.size)
sd = mc3.Gamma('sd', 1., 1., shape=locs.size)
x_obs = mc3.NormalMixture('x_obs', w, mu, sd=sd, observed=D)
trace = mc3.sample(1000, tune=4000, random_seed=SEED)[100:]
In [115]:
mc3.traceplot(trace, combined=False);
In [116]:
mc3.summary(trace)
Out[116]:
In [117]:
fig, ax = plt.subplots(3, 2, figsize=(10,12))
axs = mc3.plot_posterior(trace, point_estimate='median', ax=ax)
ax = ax.flatten()
ax[0].axvline(8000/(10000+8000), c='r', lw=3)
ax[1].axvline(1-8000/(10000+8000), c='r', lw=3)
ax[2].axvline(10, c='r', lw=3)
ax[3].axvline(3, c='r', lw=3)
ax[4].axvline(2, c='r', lw=3)
ax[5].axvline(2, c='r', lw=3)
Out[117]:
In [118]:
ppc = mc3.sample_posterior_predictive(trace, samples=205, model=model)
In [119]:
hist, bins = np.histogram(ppc['x_obs'], bins='doane', normed=True)
hist2 = []
for i in tqdm.tqdm_notebook(range(ppc['x_obs'].shape[0])):
h, _ = np.histogram(ppc['x_obs'][i], bins=bins, normed=True)
hist2.append(h)
In [120]:
hist2 = np.asarray(hist2)
hist2.shape
Out[120]:
In [121]:
perc = np.percentile(hist2, [2.5, 97.5, 50], axis=0)
perc.shape
Out[121]:
In [123]:
# plt.plot(tb.bin_edges_to_center(bins), hist2[2,:], lw=4, label='MCMC')
plt.fill_between(tb.bin_edges_to_center(bins), hist2[0,:], hist2[1,:], color='red', alpha=0.2)
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
plt.legend();
In [130]:
plt.plot(tb.bin_edges_to_center(bins), hist2.T, c='r', alpha=0.1);
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
In [93]:
np.ones_like(locs), locs.size
Out[93]:
In [94]:
SEED = 20161210
with mc3.Model() as model:
w = mc3.Dirichlet('w', np.ones(3))
mu = mc3.Uniform('mu', -10., 20., shape=3)
sd = mc3.Gamma('sd', 1., 1., shape=3)
x_obs = mc3.NormalMixture('x_obs', w, mu, sd=sd, observed=D)
trace = mc3.sample(1000, tune=4000, random_seed=SEED)[100:]
In [95]:
mc3.traceplot(trace, combined=False);
In [96]:
mc3.summary(trace)
Out[96]:
In [97]:
ppc = mc3.sample_posterior_predictive(trace, samples=205, model=model)
In [110]:
hist, bins = np.histogram(ppc['x_obs'], bins='doane', normed=True)
hist2 = []
for i in tqdm.tqdm_notebook(range(ppc['x_obs'].shape[0])):
h, _ = np.histogram(ppc['x_obs'][i], bins=bins, normed=True)
hist2.append(h)
hist2 = np.asarray(hist2)
perc = np.percentile(hist2, [2.5, 97.5, 50], axis=0)
In [111]:
hist2.shape
Out[111]:
In [112]:
plt.plot(tb.bin_edges_to_center(bins), hist2[2,:], lw=4, label='MCMC')
plt.fill_between(tb.bin_edges_to_center(bins), hist2[0,:], hist2[1,:], color='red', alpha=0.2)
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
plt.legend();
In [100]:
fig, ax = plt.subplots(5, 2, figsize=(10,12))
axs = mc3.plot_posterior(trace, point_estimate='median', ax=ax)
In [101]:
SEED = 20161210
with mc3.Model() as model:
w = mc3.Dirichlet('w', np.ones(10))
mu = mc3.Uniform('mu', -10., 20., shape=10)
sd = mc3.Gamma('sd', 1., 1., shape=10)
x_obs = mc3.NormalMixture('x_obs', w, mu, sd=sd, observed=D)
trace = mc3.sample(1000, tune=4000, random_seed=SEED)[100:]
In [102]:
ppc = mc3.sample_posterior_predictive(trace, samples=205, model=model)
In [107]:
ppc['x_obs'].shape
Out[107]:
In [108]:
hist, bins = np.histogram(ppc['x_obs'], bins='doane', normed=True)
hist2 = []
for i in tqdm.tqdm_notebook(range(ppc['x_obs'].shape[0])):
h, _ = np.histogram(ppc['x_obs'][i], bins=bins, normed=True)
hist2.append(h)
hist2 = np.asarray(hist2)
perc = np.percentile(hist2, [2.5, 97.5, 50], axis=0)
In [109]:
plt.plot(tb.bin_edges_to_center(bins), hist2[2,:], lw=4, label='MCMC')
plt.fill_between(tb.bin_edges_to_center(bins), hist2[0,:], hist2[1,:], color='red', alpha=0.2)
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
plt.legend();
In [106]:
mc3.summary(trace)
Out[106]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
with mc3.Model() as model:
w = mc3.Dirichlet('w', np.ones_like(locs))
mu1 = mc3.Uniform('mu1', -10., 20.)
mu2 = mc3.Uniform('mu2', mu1, 20.)
sd = mc3.Gamma('sd', 1., 1., shape=locs.size)
x_obs = mc3.NormalMixture('x_obs', w, [mu1, mu2], sd=sd, observed=D)
start = mc3.find_MAP()
trace = mc3.sample(5000, n_init=10000, tune=1000, random_seed=SEED, start=start)[1000:]
In [ ]:
mc3.traceplot(trace, combined=True);
mc3.summary(trace)
In [ ]:
ppc = mc3.sample_ppc(trace, samples=5000, model=model)
hist, bins = np.histogram(ppc['x_obs'], bins='doane', normed=True)
plt.plot(tb.bin_edges_to_center(bins), hist, lw=4, label='MCMC')
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
plt.legend()
In [ ]:
D.shape
In [ ]:
with mc3.Model() as model:
w = mc3.Dirichlet('w', np.ones_like(locs))
mu = mc3.Uniform('mu', -10., 20., shape=locs.size)
sd = mc3.Gamma('sd', 1., 1., shape=locs.size)
x_obs = mc3.Mixture('x_obs', w=w, comp_dists=mc3.Normal.dist(mu=mu, sd=sd), observed=D)
start = mc3.find_MAP()
trace = mc3.sample(5000, n_init=10000, tune=1000, random_seed=SEED, start=start)[1000:]
In [ ]:
def test_mixture_list_of_poissons(self):
with Model() as model:
w = Dirichlet('w', np.ones_like(self.pois_w))
mu = Gamma('mu', 1., 1., shape=self.pois_w.size)
x_obs = Mixture('x_obs', w,
[Poisson.dist(mu[0]), Poisson.dist(mu[1])],
observed=self.pois_x)
step = Metropolis()
trace = sample(5000, step, random_seed=self.random_seed, progressbar=False)
assert_allclose(np.sort(trace['w'].mean(axis=0)),
np.sort(self.pois_w),
rtol=0.1, atol=0.1)
assert_allclose(np.sort(trace['mu'].mean(axis=0)),
np.sort(self.pois_mu),
rtol=0.1, atol=0.1)
In [ ]:
ppc = mc3.sample_ppc(trace, samples=5000, model=model)
hist, bins = np.histogram(ppc['x_obs'], bins='doane', normed=True)
plt.plot(tb.bin_edges_to_center(bins), hist, lw=4, label='MCMC')
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
plt.legend()
In [ ]:
mc3.traceplot(trace, combined=True);
mc3.summary(trace)
In [ ]:
In [ ]:
with mc3.Model() as model:
w = mc3.Dirichlet('w', np.ones_like(locs))
mu1 = mc3.Uniform('mu1', -10., 20.,)
mu2 = mc3.Uniform('mu2', mu1, 20., )
sd = mc3.Gamma('sd', 1., 1., shape=locs.size)
x_obs = mc3.Mixture('x_obs', w=w, comp_dists=mc3.Normal.dist(mu=[mu1,mu2], sd=sd, shape=locs.size),
observed=D)
start = mc3.find_MAP()
trace = mc3.sample(5000, n_init=10000, tune=1000, random_seed=SEED, start=start)[1000:]
In [ ]:
mc3.traceplot(trace, combined=True);
mc3.summary(trace)
In [ ]:
ppc = mc3.sample_ppc(trace, samples=5000, model=model)
hist, bins = np.histogram(ppc['x_obs'], bins='doane', normed=True)
plt.plot(tb.bin_edges_to_center(bins), hist, lw=4, label='MCMC')
plt.hist(D, bins=bins, histtype='step', lw=3, normed=True, label='Data');
plt.legend()
In [ ]:
In [ ]:
In [ ]: